Mars2020_DataAnalysis

IMPORT
Librairies
Data
ABUNDANCES
NORMALIZE DATA
COMPUTE DISTANCES
- 1. UniFrac, Aitchison, Bray-Curtis and Canberra
- 2. Plot in 3D
HIERARCHICAL CLUSTERING
REPRODUCE PLOTS FROM PAPER

IMPORT

Librairies

library(phyloseq) # for phyloseq object
library(ggplot2)
library(cowplot)
library(plyr)
library(dplyr)
library("plotly") # plot 3D
library("microbiome") # for centered log-ratio
library("coda") # Aitchison distance
library("coda.base") # Aitchison distance
library("vegan") # NMDS
library(pheatmap) # for heatmap

Data

# Set path
path <- "~/Projects/IBS_Meta-analysis_16S"

# Import phyloseq object
physeq.mars <- readRDS(file.path(path, "phyloseq-objects/physeq_mars.rds"))
physeq.mars <- prune_samples(sample_sums(physeq.mars)>=500, physeq.mars)

# Sanity check
physeq.mars

## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 1561 taxa and 69 samples ]
## sample_data() Sample Data:       [ 69 samples by 12 sample variables ]
## tax_table()   Taxonomy Table:    [ 1561 taxa by 7 taxonomic ranks ]

Phylogenetic tree was computed with the package phangorn, and the script was run on a cluster. Let’s check we have correctly generated a phylogenetic tree.

# Look at the tree
plot_tree(physeq.mars, color = "host_disease", ladderize="left")

ABUNDANCES

1. Absolute abundances

# Plot Phylum
plot_bar(physeq.mars, fill = "Phylum") + facet_wrap("host_disease", scales="free") +
  theme(axis.text.x = element_text(size = 8))+
  labs(x = "Samples", y = "Absolute abundance", title = "Mars dataset (2020)")

# Plot Class
plot_bar(physeq.mars, fill = "Class")+ facet_wrap("host_disease", scales="free") +
  theme(axis.text.x = element_text(size = 8))+
  labs(x = "Samples", y = "Absolute abundance", title = "Mars dataset (2020)")

Sequencing depth characteristics of the Mars dataset:
- minimum of 1140 total count per sample
- median: 3.255610^{4} total count per sample
- maximum of 6.360210^{4} total count per sample

2. Relative abundances

# Agglomerate to phylum & class levels
phylum.table <- physeq.mars %>%
  tax_glom(taxrank = "Phylum") %>%                     # agglomerate at phylum level
  transform_sample_counts(function(x) {x/sum(x)} ) %>% # Transform to rel. abundance
  psmelt()                                             # Melt to long format

class.table <- physeq.mars %>%
  tax_glom(taxrank = "Class") %>%
  transform_sample_counts(function(x) {x/sum(x)} ) %>%
  psmelt()


# Plot relative abundances
ggplot(phylum.table, aes(x = Sample, y = Abundance, fill = Phylum))+
  facet_wrap(~ host_disease, scales = "free") + # scales = "free" removes empty lines
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(size = 8, angle = -90))+
  labs(x = "Samples", y = "Relative abundance", title = "Mars dataset (2020)")

ggplot(class.table, aes(x = Sample, y = Abundance, fill = Class))+
  facet_wrap(~ host_disease, scales = "free") + # scales = "free" removes empty lines
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(size = 8, angle = -90))+
  labs(x = "Samples", y = "Relative abundance", title = "Mars dataset (2020)")

3. Firmicutes/Bacteroidota ratio

# Extract abundance of only Bacteroidota and Firmicutes
bacter <- phylum.table %>%
  filter(Phylum == "Bacteroidota") %>%
  select(c('Sample', 'Abundance', 'host_disease', 'Phylum', 'host_subtype', 'Collection')) %>%
  arrange(Sample)

firmi <- phylum.table %>%
  filter(Phylum == "Firmicutes") %>%
  select(c('Sample', 'Abundance', 'host_disease', 'Phylum', 'host_subtype', 'Collection')) %>%
  arrange(Sample)

# Calculate log2 ratio Firmicutes/Bacteroidota
ratio.FB <- data.frame('Sample' = bacter$Sample,
                       'host_disease' = bacter$host_disease,
                       'host_subtype' = bacter$host_subtype,
                       'Collection' = bacter$Collection,
                       'Bacteroidota' = bacter$Abundance,
                       'Firmicutes' = firmi$Abundance)
ratio.FB$logRatioFB <- log2(ratio.FB$Firmicutes / ratio.FB$Bacteroidota)

# Plot log2 ratio Firmicutes/Bacteroidota
ggplot(ratio.FB, aes(x = host_disease, y = logRatioFB))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(width=0.1)+
  labs(x = "",  y = 'Log2(Firmicutes/Bacteroidota)', title = "Firmicutes:Bacteroidota ratio")+
  theme_cowplot()

# Plot by IBS subtype
ggplot(ratio.FB, aes(x = host_subtype, y = logRatioFB))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(width=0.1)+
  labs(x = "",  y = 'Log2(Firmicutes/Bacteroidota)', title = "Firmicutes:Bacteroidota ratio")+
  theme_cowplot()

# Plot by IBS subtype with collection time 1 or 2
ggplot(ratio.FB, aes(x = host_subtype, y = logRatioFB))+
  geom_boxplot(outlier.shape = NA)+
  geom_jitter(width=0.1)+
  facet_wrap(~Collection)+
  labs(x = "",  y = 'Log2(Firmicutes/Bacteroidota)', title = "Firmicutes:Bacteroidota ratio")+
  theme_cowplot()+
  theme(axis.text.x = element_text(angle=90))

NORMALIZE DATA

# Sanity check no sample with less than 500 total count
table(sample_sums(physeq.mars)<500) # all FALSE

#____________________________________________________________________
# PHYLOSEQ OBJECT WITH NON-ZERO COMPOSITIONS
physeq.NZcomp <- physeq.mars
otu_table(physeq.NZcomp)[otu_table(physeq.NZcomp) == 0] <- 0.5 # pseudocounts

# Sanity check that 0 values have been replaced
# otu_table(physeq.mars)[1:5,1:5]
# otu_table(physeq.NZcomp)[1:5,1:5]

# transform into compositions
physeq.NZcomp <- transform_sample_counts(physeq.NZcomp, function(x) x / sum(x) )
table(rowSums(otu_table(physeq.NZcomp))) # check if there is any row not summing to 1

# Save object
saveRDS(physeq.NZcomp, file.path(path, "data/analysis-individual/Mars-2020/02_EDA-Mars/physeq_NZcomp.rds"))

#____________________________________________________________________
# PHYLOSEQ OBJECT WITH RELATIVE COUNT (BETWEEN 0 AND 1)
physeq.rel <- physeq.mars
physeq.rel <- transform_sample_counts(physeq.rel, function(x) x / sum(x) ) # divide each count by the total number of counts (per sample)

# check the counts are all relative
# otu_table(physeq.mars)[1:5, 1:5]
# otu_table(physeq.rel)[1:5, 1:5]

# sanity check
table(rowSums(otu_table(physeq.rel))) # check if there is any row not summing to 1

# save the physeq.rel object
saveRDS(physeq.rel, file.path(path, "data/analysis-individual/Mars-2020/02_EDA-Mars/physeq_relative.rds"))

#____________________________________________________________________
# PHYLOSEQ OBJECT WITH COMMON-SCALE NORMALIZATION
physeq.CSN <- physeq.mars
physeq.CSN <- transform_sample_counts(physeq.CSN, function(x) (x*min(sample_sums(physeq.CSN))) / sum(x) )

# sanity check
table(rowSums(otu_table(physeq.CSN))) # check that all rows are summing to the same total

# save the physeq.CSN object
saveRDS(physeq.CSN, file.path(path, "data/analysis-individual/Mars-2020/02_EDA-Mars/physeq_CSN.rds"))


#____________________________________________________________________
# PHYLOSEQ OBJECT WITH CENTERED LOG RATIO COUNT
physeq.clr <- physeq.mars
physeq.clr <- microbiome::transform(physeq.mars, "clr") # the function adds pseudocounts itself

# Compare the otu tables in the original phyloseq object and the new one after CLR transformation
otu_table(physeq.mars)[1:5, 1:5] # should contain absolute counts
otu_table(physeq.clr)[1:5, 1:5] # should all be relative

# save the physeq.rel object
saveRDS(physeq.clr, file.path(path, "data/analysis-individual/Mars-2020/02_EDA-Mars/physeq_clr.rds"))

COMPUTE DISTANCES

1. UniFrac, Aitchison, Bray-Curtis and Canberra

First, let’s look at these four distances of interest.

#____________________________________________________________________________________
# Measure distances
getDistances <- function(){
  # set.seed(123) # for unifrac, need to set a seed
  # glom.UniF <- UniFrac(physeq.rel, weighted=TRUE, normalized=TRUE) # weighted unifrac
  glom.ait <- phyloseq::distance(physeq.clr, method = 'euclidean') # aitchison
  glom.bray <- phyloseq::distance(physeq.CSN, method = "bray") # bray-curtis
  glom.can <- phyloseq::distance(physeq.NZcomp, method = "canberra") # canberra
  dist.list <- list(#"UniF" = glom.UniF,
    "Ait" = glom.ait, "Canb" = glom.can, "Bray" = glom.bray)
  
  return(dist.list)
}


#____________________________________________________________________________________
# Plot in 2D the distances
plotDistances2D <- function(dlist, ordination="MDS"){
  plist <- NULL
  # plist <- vector("list", 4)
  # names(plist) <- c("Weighted Unifrac", "Aitchison", "Bray-Curtis", "Canberra")
  plist <- vector("list", 3)
  names(plist) <- c("Aitchison", "Bray-Curtis", "Canberra")
  
  # print("Unifrac")
  # # Weighted UniFrac
  # set.seed(123)
  # iMDS.UniF <- ordinate(physeq.rel, ordination, distance=dlist$UniF)
  # plist[[1]] <- plot_ordination(physeq.rel, iMDS.UniF, color="host_disease")
  
  print("Aitchison")
  # Aitchison
  set.seed(123)
  iMDS.Ait <- ordinate(physeq.clr, ordination, distance=dlist$Ait)
  plist[[1]] <- plot_ordination(physeq.clr, iMDS.Ait, color="host_disease")
  
  print("Bray")
  # Bray-Curtis
  set.seed(123)
  iMDS.Bray <- ordinate(physeq.CSN, ordination, distance=dlist$Bray)
  plist[[2]] <- plot_ordination(physeq.CSN, iMDS.Bray, color="host_disease")
  
  print("Canberra")
  # Canberra
  set.seed(123)
  iMDS.Can <- ordinate(physeq.NZcomp, ordination, distance=dlist$Can)
  plist[[3]] <- plot_ordination(physeq.NZcomp, iMDS.Can, color="host_disease")
  
  # Creating a dataframe to plot everything
  plot.df = ldply(plist, function(x) x$data)
  names(plot.df)[1] <- "distance"
  
  return(plot.df)
}

Now let’s plot!

# Get the distances & the plot data
dist.mars <- getDistances()
plot.df <- plotDistances2D(dist.mars)

## [1] "Aitchison"
## [1] "Bray"
## [1] "Canberra"

# Plot
ggplot(plot.df, aes(Axis.1, Axis.2, color=host_disease))+
  geom_point(size=6, alpha=0.5)  + scale_color_manual(values = c('blue', 'red'))+
  facet_wrap(distance~., scales='free', nrow=1)+
  theme_bw()+
  theme(strip.text.x = element_text(size=20))+
  labs(color="Disease")

# ggsave(file.path(path.plots, "distances4_MDS.jpg"), height = 4, width = 15)
# ggsave(file.path(path.plots, "distances4_MDS.jpg"), height = 4, width = 10)

2. Plot in 3D

For better visualization, we will also take a glance at reduction to 3D.

#____________________________________________________________________________________
# Plot 3D ordination
plotDistances3D <- function(d, name_dist){
  
  # Reset parameters
  mds.3D <- NULL
  xyz <- NULL
  fig.3D <- NULL
  
  # Reduce distance matrix to 3 dimensions
  set.seed(123)
  mds.3D <- metaMDS(d, method="MDS", k=3, trace = 0)
  xyz <- scores(mds.3D, display="sites") # pull out the (x,y,z) coordinates
  
  # Plot
  fig.3D <- plot_ly(x=xyz[,1], y=xyz[,2], z=xyz[,3], type="scatter3d", mode="markers",
                    color=sample_data(physeq.mars)$host_disease, colors = c("blue", "red"))%>%
    layout(title = paste('MDS in 3D with', name_dist, 'distance', sep = ' '))
  
  return(fig.3D)
}